#pragma once

#include "InterpCondition.h"
#include "Solver.h"
#include "Vec.h"

#define MAX(a, b) ((a > b) ? a : b)
#define MIN(a, b) ((a > b) ? b : a)

// 关于Vec的宏
#define Vec1 Vec<double, 1>
#define Vec2 Vec<double, 2>
#define VEC(Dim) Vec<double, (Dim)>

namespace SplineSp
{
    using namespace std;

    // 差分表类型
    using Table = vector<vector<double>>;

    enum class SplineType
    {
        ppForm,   // 一般样条
        cardinalB // 基数 B 样条
    };

    enum class BCType
    {
        complete, // 完全样条
        notAknot, // 无结样条
        periodic  // 周期样条
    };

    // 多项式
    template <int Order, class CoefType>
    class Polynomial
    {
    private:
        CoefType coefficients[Order + 1]; //系数向量

    public:
        Polynomial()
        {
            for (int i = 0; i < Order + 1; i++)
            {
                coefficients[i] = 0;
            }
        }
        Polynomial(const CoefType arr[])
        {
            for (int i = 0; i < Order + 1; i++)
            {
                coefficients[i] = arr[i];
            }
        }
        Polynomial(std::initializer_list<CoefType> l)
        {
            auto it = l.begin();
            for (int i = 0; i < Order + 1; i++)
            {
                coefficients[i] = *it++;
            }
        }
        Polynomial(const Polynomial &poly)
        {
            for (int i = 0; i < Order + 1; i++)
            {
                coefficients[i] = poly[i];
            }
        }

        // 多项式赋值，注意考虑不同次系数向量长度不同的问题，需要用 Polynomial() 初始化
        template <int Order2, class CoefType2>
        Polynomial(const Polynomial<Order2, CoefType2> &poly) : Polynomial()
        {
            for (int i = 0; i < MIN(Order + 1, Order2 + 1); i++)
            {
                coefficients[i] = poly[i];
            }
        }

// 多项式加法
#define POLY_BINARY_OP(OpName, Op)                              \
    template <int Order2>                                       \
    auto OpName(const Polynomial<Order2, CoefType> &poly) const \
    {                                                           \
        int index = 0;                                          \
        int this_size = order() + 1;                            \
        int that_size = poly.order() + 1;                       \
        CoefType coef[MAX(this_size, that_size)];               \
        while (index < MIN(this_size, that_size))               \
        {                                                       \
            coef[index] = coefficients[index] Op poly[index];   \
            index++;                                            \
        }                                                       \
        while (index < that_size)                               \
        {                                                       \
            int first = 0;                                      \
            first = first Op 1;                                 \
            coef[index] = poly[index] * first;                  \
            index++;                                            \
        }                                                       \
        return Polynomial<MAX(Order, Order2), CoefType>(coef);  \
    }

        POLY_BINARY_OP(operator+, +)
        POLY_BINARY_OP(operator-, -)
#undef POLY_BINARY_OP

// 多项式乘系数
#define POLY_NUM_OP(OpName, Op)             \
    Polynomial OpName(CoefType k) const     \
    {                                       \
        CoefType coef[Order + 1];           \
        for (int i = 0; i < Order + 1; i++) \
        {                                   \
            coef[i] = coefficients[i] Op k; \
        }                                   \
        return Polynomial(coef);            \
    }

        POLY_NUM_OP(operator*, *)
        POLY_NUM_OP(operator/, /)
#undef POLY_NUM_OP

        // 多项式相乘
        template <int Order2>
        auto operator*(const Polynomial<Order2, CoefType> &poly)
        {
            const int m = Order + 1;
            const int n = Order2 + 1;
            CoefType coef[m + n];
            // 注意清空数组
            for (int i = 0; i < m + n; i++)
            {
                coef[i] = 0;
            }
            for (int i = 0; i < m; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    coef[i + j] = coef[i + j] + coefficients[i] * poly[j];
                }
            }
            return Polynomial<Order + Order2, CoefType>(coef);
        }

        CoefType &operator[](int i) { return coefficients[i]; }
        CoefType operator[](int i) const { return coefficients[i]; }

        // 向右平移多项式 dx
        Polynomial<Order, CoefType> operator>>(double dx)
        {
            Polynomial<Order, CoefType> origin, poly;
            origin[0] = -1 * dx, origin[1] = 1;
            for (int i = 0; i < Order + 1; i++)
            {
                Polynomial<Order, CoefType> tmp;
                tmp[0] = 1;
                for (int j = 0; j < i; j++)
                {
                    tmp = tmp * origin;
                }
                poly = poly + tmp * coefficients[i];
            }
            return poly;
        }

        friend ostream &operator<<(ostream &out, const Polynomial &poly)
        {
            out << poly.coefficients[0];

            for (int i = 1; i < (int)poly.order() + 1; i++)
            {
                out << " + " << poly[i] << "x^" << i;
            }
            out << endl;
            return out;
        }

        // 求值
        CoefType operator()(double x)
        {
            CoefType y = 0;
            double x_i = 1;
            for (int i = 0; i < Order + 1; i++)
            {
                y = y + coefficients[i] * x_i;
                x_i = x_i * x;
            }
            return y;
        }

        // 求函数的幂
        template <int T, bool stop>
        Polynomial<Order * T, CoefType> pow()
        {
            Polynomial<Order * T, CoefType> poly(coefficients);
            Polynomial<Order * T, CoefType> origin = poly;
            for (int i = 1; i < T; i++)
            {
                Polynomial<Order * T, CoefType> tmp = poly * origin;
                poly = tmp;
            }
            return poly;
        }

        // 求导
        Polynomial<Order - 1, CoefType> derivation()
        {
            CoefType _coefficients[Order];
            for (int i = 1; i < Order + 1; i++)
            {
                _coefficients[i - 1] = coefficients[i] * i;
            }
            return Polynomial<Order - 1, CoefType>(_coefficients);
        }

        // 绘制图像
        void plot(double left, double right, char color)
        {
            //调整输出精度
            cout << setprecision(16);
            cout << "x = linspace(" << left << "," << right << ");" << endl;
            cout << "plot(x," << coefficients[0];

            for (int i = 1; i < Order + 1; i++)
            {
                cout << " + " << coefficients[i] << "*x.^" << i;
            }
            cout << ",\'" << color << "\');" << endl;
            cout << setprecision(6);
        }

        // 绘制二维图像
        void plot2(double left, double right, char color)
        {
            cout << setprecision(16);
            cout << "t = linspace(" << left << "," << right << ");" << endl;
            cout << "plot(" << coefficients[0][0];

            for (int i = 1; i < Order + 1; i++)
            {
                cout << " + " << coefficients[i][0] << "*t.^" << i;
            }

            cout << ", " << coefficients[0][1];

            for (int i = 1; i < Order + 1; i++)
            {
                cout << " + " << coefficients[i][1] << "*t.^" << i;
            }
            cout << ",\'" << color << "\');" << endl;
            cout << setprecision(6);
        }

        const int order() const { return Order; } //获取多项式的阶
    };

    /**
     * @date   2021/11/24 11:43:50
     *
     * @arg    const Polynomial<Order, VEC(Dim1)> &poly1, const Polynomial<Order, VEC(Dim2)> &poly2
     * @return Polynomial<Order, VEC(Dim1 + Dim2)>
     *
     */

    //将两个同阶多项式合并为高维多项式
    template <int Order, int Dim1, int Dim2>
    auto overlying(const Polynomial<Order, VEC(Dim1)> &poly1, const Polynomial<Order, VEC(Dim2)> &poly2)
    {
        Vec<double, Dim1 + Dim2> coef[Order + 1];
        for (int i = 0; i < Order + 1; i++)
        {
            for (int j = 0; j < Dim1; j++)
            {
                coef[i][j] = poly1[i][j];
            }
            for (int j = Dim1; j < Dim1 + Dim2; j++)
            {
                coef[i][j] = poly2[i][j - Dim1];
            }
        }
        return Polynomial<Order, VEC(Dim1 + Dim2)>(coef);
    }

    // ================================================================================

    template <int Dim, int Order, SplineType t>
    class Spline
    {
    private:
        using CoefType = Vec<double, Dim>;              // 规定系数类型
        vector<double> interval;                        // 非零区间，长度比 splinePoly 多 1
        vector<Polynomial<Order, CoefType>> splinePoly; // 分段多项式

    public:
        Spline(int len, const double arr[]) : splinePoly(vector<Polynomial<Order, CoefType>>(len - 1))
        {
            for (int i = 0; i < len; i++)
            {
                interval.push_back(arr[i]);
            }
        }
        Spline(std::initializer_list<double> l) : splinePoly(vector<Polynomial<Order, CoefType>>(l.size() - 1))
        {
            auto it = l.begin();
            for (int i = 0; i < l.size(); i++)
            {
                interval.push_back(*it++);
            }
        }
        Spline(vector<double> l) : splinePoly(vector<Polynomial<Order, CoefType>>(l.size() - 1))
        {
            interval = l;
        }
        Spline(const Spline &spline)
        {
            splinePoly = spline.splinePoly;
            interval = spline.interval;
        }

        // 导数
        Spline<Dim, Order - 1, t> derivation()
        {
            Spline<Dim, Order - 1, t> res(interval);
            for (int i = 0; i < size(); i++)
            {
                res[i] = splinePoly[i].derivation();
            }
            return res;
        }

        // 计算值
        CoefType operator()(double x) { return getPoly(x)(x); }

        // 取分段多项式
        Polynomial<Order, CoefType> operator[](int index) const { return (index < size() ? splinePoly[index] : splinePoly[0]); }
        Polynomial<Order, CoefType> &operator[](int index) { return (index < size() ? splinePoly[index] : splinePoly[0]); }

        // 根据点所在区间取多项式
        // 注意考虑左闭右开区间
        Polynomial<Order, CoefType> getPoly(double x) const
        {
            for (int i = 0; i < size(); i++)
            {
                // 按左开右闭区间
                if (x >= interval[i] && x < interval[i + 1])
                {
                    return splinePoly[i];
                }
            }
            return Polynomial<Order, CoefType>();
        }

// 基数B样条加法：注意未初始化的样条
// 注意调整区间和多项式向量的长度
#define SPLINE_BINARY_OP(OpName, Op)                                                \
    template <int Order2>                                                           \
    auto OpName(const Spline<1, Order2, SplineType::cardinalB> &spline2)            \
    {                                                                               \
        int l = MIN(left(), spline2.left());                                        \
        int r = MAX(right(), spline2.right());                                      \
        double inter[r + 1 - l];                                                    \
        for (int i = l; i < r + 1; i++)                                             \
        {                                                                           \
            inter[i - l] = i;                                                       \
        }                                                                           \
        Spline<1, MAX(Order, Order2), SplineType::cardinalB> res(r + 1 - l, inter); \
        for (int i = l; i < r; i++)                                                 \
        {                                                                           \
            res[i - l] = getPoly(i) Op spline2.getPoly(i);                          \
        }                                                                           \
        return res;                                                                 \
    }

        SPLINE_BINARY_OP(operator+, +)
        SPLINE_BINARY_OP(operator-, -)
#undef SPLINE_BINARY_OP

// 样条乘法
#define SPLINE_NUM_OP(OpName, Op)                        \
    Spline<Dim, Order, t> OpName(const CoefType k) const \
    {                                                    \
        Spline<Dim, Order, t> res = *this;               \
        for (int i = 0; i < size(); i++)                 \
        {                                                \
            res[i] = res[i] Op k;                        \
        }                                                \
        return res;                                      \
    }

        SPLINE_NUM_OP(operator*, *)
        SPLINE_NUM_OP(operator/, /)
#undef SPLINE_NUM_OP

        // 多项式乘样条
        template <int Order2>
        // 传参时如果不改变参数值，尽可能把参数写成const，防止出现传入常量时编译产生的问题
        Spline<Dim, Order + Order2, t> operator*(const Polynomial<Order2, CoefType> &poly)
        {
            Spline<Dim, Order + Order2, t> res(interval);
            for (int i = 0; i < size(); i++)
            {
                res[i] = splinePoly[i] * poly;
            }
            return res;
        }

        // 样条赋值
        template <int Order2, SplineType t2>
        void operator=(const Spline<Dim, Order2, t2> &spline)
        {
            interval = spline.getInterval();
            splinePoly.clear();
            for (int i = 0; i < spline.size(); i++)
            {
                splinePoly.push_back(spline[i]);
            }
        }

        // 右移样条
        Spline<Dim, Order, t> operator>>(double k)
        {
            Spline<Dim, Order, t> res = *this;
            res.interval[0] += k;
            for (int i = 0; i < size(); i++)
            {
                res[i] = res[i] >> k;
                res.interval[i + 1] += k;
            }
            return res;
        }

        // 标准输出
        friend ostream &operator<<(ostream &out, const Spline &spline)
        {
            for (int i = 0; i < spline.size(); i++)
            {
                out << "(" << spline.interval[i] << "," << spline.interval[i + 1] << ") "
                    << "s" << i + 1 << "(x) = " << spline[i];
            }
            return out;
        }

        //绘制图像
        void plot(char color)
        {
            for (int i = 0; i < splinePoly.size(); i++)
            {
                splinePoly[i].plot(interval[i], interval[i + 1], color);
                cout << "hold on;" << endl;
            }
        }

        void plot2(char color)
        {
            for (int i = 0; i < splinePoly.size(); i++)
            {
                splinePoly[i].plot2(interval[i], interval[i + 1], color);
                cout << "hold on;" << endl;
            }
        }

        // 获取和修改成员变量
        int dim() const { return Dim; }                         // 系数的维数
        int order() const { return Order; }                     // 多项式阶数
        int size() const { return splinePoly.size(); }          // 分段数
        double left() const { return interval.front(); }        // 区间左端点
        double right() const { return interval.back(); }        // 区间右端点
        vector<double> getInterval() const { return interval; } // 获取区间

        // =========================================================================

        // 标准样条插值
        template <int Ord, SplineType _t>
        friend Spline<1, Ord, _t> interpolate(const InterpC &, BCType);

        // 计算部分差分表，宽度为 Order
        void calculateTable(int width, InterpC &interpC, Table &);

        // 逐段计算样条，同时通过每段多项式添加新的边值条件
        Polynomial<Order, double> pieceCalculate(InterpC &interpC, vector<double> &M, int index);

        // 基数B样条
        template <int Ord>
        friend Spline<1, Ord, SplineType::cardinalB> BSpline(int i);
    };

    template <int Dim, int Order, SplineType t>
    void Spline<Dim, Order, t>::calculateTable(int width, InterpC &interpC, Table &table)
    {
        vector<double> cond;
        // 计算重数地将插值点推入
        for (int i = 0; i < interval.size(); i++)
        {
            for (int j = 0; j < interpC.size(i); j++)
            {
                cond.push_back(interval[i]);
            }
        }
        // 对每个插值点计算
        for (int i = 0; i < cond.size(); i++)
        {
            // 如果是第一个插值条件，就直接推入
            table.push_back({interpC(cond[i])});
            // 注意边界条件取最小值
            for (int j = 1; j < MIN(width, i + 1); j++)
            {
                double left = cond[i - j], right = cond[i];
                // 如果左右不相等，就计算两点差分
                if (left != right)
                {
                    table[i].push_back((table[i][j - 1] - table[i - 1][j - 1]) / (right - left));
                    continue;
                }
                // 否则推入该点的 factor 阶导数
                int factor = 1;
                for (int k = 1; k <= j; factor *= k, k++)
                    ;
                table[i].push_back(1.0 / factor * interpC(cond[i], factor));
            }
        }
    }

    // 逐段计算多项式
    template <int Dim, int Order, SplineType t>
    Polynomial<Order, double> Spline<Dim, Order, t>::pieceCalculate(InterpC &interpC, vector<double> &M, int index)
    {
        Polynomial<Order, double> res;
        // 泰勒展开的位置
        double x = interpC[index];
        Polynomial<0, double> res_0({interpC(x)});
        Polynomial<1, double> res_1({-1 * x, 1});
        Polynomial<2, double> res_2 = res_1 * res_1;
        Polynomial<3, double> res_3 = res_1 * res_2;

        // 利用泰勒展开式
        // s_i(x) = f_i + f[x_i, x_i](x - x_i) + M_i / 2 * (x - x_i)^2 + (M_{i+1} - M_i) / 6(x_{i+1} - x_i) * (x - x_i)^3
        res = res_0 + res_1 * interpC(x, 1) + res_2 * M[index] / 2 +
              res_3 * (M[index + 1] - M[index]) / (6 * (interpC[index + 1] - interpC[index]));

        // 添加新的边值条件
        double m_1 = res.derivation()(interpC[index + 1]);
        interpC.sites(index + 1).push_back(m_1);

        return res;
    }

    // 计算基数B样条
    template <int Ord>
    Spline<1, Ord, SplineType::cardinalB> BSpline(int i)
    {
        Spline<1, Ord, SplineType::cardinalB> res({i - 1.0, i * 1.0});
        res[0] = Polynomial<0, Vec1>({1});

        // 利用递推公式和平移性质
        for (int j = 0; j < Ord; j++)
        {
            // 注意每次阶数 Order = j + 1
            res = res * Polynomial<1, Vec1>({(1.0 - i) / (j + 1), 1.0 / (j + 1)}) +
                  (res >> 1) * Polynomial<1, Vec1>({1.0 * (j + 1 + i) / (j + 1), -1.0 / (j + 1)});
        }
        return res;
    }

    // 标准三类样条插值和二次、三次B样条
    template <int Ord, SplineType t>
    Spline<1, Ord, t> interpolate(InterpC &interpC, BCType type)
    {
        Spline<1, Ord, t> res(interpC.points); // 将插值条件作为区间初始化
        Table table;                           // 差分表
        int N = interpC.size();                // 插值点的个数

        // 根据不同类型的插值计算
        switch (t)
        {
        case SplineType::ppForm:
        {
            res.calculateTable(Ord, interpC, table); // 计算差分表

            // 线性系统
            vector<vector<double>> A(N, vector<double>(N));
            vector<double> B(N), lam(N - 2), mu(N - 2);

            // 计算系数
            for (int i = 0; i < N - 2; i++)
            {
                mu[i] = (interpC[i + 1] - interpC[i]) / (interpC[i + 2] - interpC[i]);
                lam[i] = (interpC[i + 2] - interpC[i + 1]) / (interpC[i + 2] - interpC[i]);
            }
            // 系数矩阵
            for (int i = 1; i < N - 1; i++)
            {
                A[i][i - 1] = mu[i - 1];
                A[i][i + 1] = lam[i - 1];
                A[i][i] = 2;
            }

            //完全样条
            if (type == BCType::complete)
            {
                // 构造线性方程组
                A[0][0] = 2, A[0][1] = 1;
                A[N - 1][N - 2] = 1, A[N - 1][N - 1] = 2;

                for (int i = 0; i < N; i++)
                {
                    B[i] = 6 * table[i + 2][Ord - 1];
                }
            }
            // 无结样条
            if (type == BCType::notAknot)
            {
                // 利用三阶导相等
                A[0][0] = interpC[2] - interpC[1], A[0][1] = interpC[0] - interpC[2], A[0][2] = interpC[1] - interpC[0];
                A[N - 1][N - 3] = interpC[N - 1] - interpC[N - 2], A[N - 1][N - 2] = interpC[N - 3] - interpC[N - 1], A[N - 1][N - 1] = interpC[N - 2] - interpC[N - 3];

                for (int i = 1; i < N - 1; i++)
                {
                    B[i] = 6 * table[i + 1][Ord - 1];
                }
                B[0] = B[N - 1] = 0;
            }
            // 周期样条
            if (type == BCType::periodic)
            {
                // 构造线性方程组
                A[0][0] = 1, A[0][N - 1] = -1;
                A[N - 1][0] = 2 * (interpC[1] - interpC[0]), A[N - 1][1] = interpC[1] - interpC[0];
                A[N - 1][N - 2] = interpC[N - 1] - interpC[N - 2], A[N - 1][N - 1] = 2 * (interpC[N - 1] - interpC[N - 2]);

                for (int i = 1; i < N - 1; i++)
                {
                    B[i] = 6 * table[i + 1][Ord - 1];
                }
                // 6(f[x_1, x_2] - f[x_N - x_{N-1}])
                B[0] = 0, B[N - 1] = 6 * (interpC(interpC[1]) - interpC(interpC[0])) / (interpC[1] - interpC[0]) -
                                     (interpC(interpC[N - 1]) - interpC(interpC[N - 2])) / (interpC[N - 1] - interpC[N - 2]);
            }

            // 求解二阶导的向量
            vector<double> M = Solver::solve(A, B);

            // 无结样条/周期样条
            if (type == BCType::notAknot || type == BCType::periodic)
            {
                // 根据二阶导向量补充一阶导的信息
                // f[x_1, x_1] = f[x_1, x_2] - (2M_1 + M_2) / 6 * (x_2 - x_1)
                double fPrime = table[1][1] - (2 * M[0] + M[1]) / 6 * (interpC[1] - interpC[0]);
                interpC.sites(0).push_back(fPrime);
            }

            // 利用导数条件推导多项式
            for (int i = 0; i < N - 1; i++)
            {
                // 将 double 转换为 Vec1
                res[i] = Polynomial<Ord, Vec1>(res.pieceCalculate(interpC, M, i));
            }

            break;
        }
        case SplineType::cardinalB:
        {
            // 三次 B 样条，在整数点插值
            if (Ord == 3 && type == BCType::complete)
            {
                vector<double> B(N);
                vector<vector<double>> M(N, vector<double>(N));

                // 系数多项式
                for (int i = 1; i < N - 1; i++)
                {
                    M[i][i - 1] = M[i][i + 1] = 1;
                    M[i][i] = 4;
                }
                M[0][0] = M[N - 1][N - 1] = 4;
                M[0][1] = M[N - 1][N - 2] = 2;

                // 目标向量
                for (int i = 1; i < N - 1; i++)
                {
                    B[i] = 6 * interpC(interpC[i]);
                }
                B[0] = 6 * interpC(interpC[0]) + 2 * interpC(interpC[0], 1);
                B[N - 1] = 6 * interpC(interpC[N - 1]) - 2 * interpC(interpC[N - 1], 1);

                // 求解获得样条系数
                vector<double> a(N + 2);
                vector<double> tmp_a = Solver::solve(M, B);

                for (int i = 0; i < N; i++)
                {
                    a[i + 1] = tmp_a[i];
                }
                a[0] = a[2] - 2 * interpC(interpC[0], 1);
                a[N + 1] = a[N - 1] + 2 * interpC(interpC[N - 1], 1);

                // 累计样条函数
                res = BSpline<3>(-1) * a[0];
                for (int i = 0; i < N + 1; i++)
                {
                    res = res + BSpline<3>(i) * a[i + 1];
                }
                // 根据左端点的位置平移
                res = (res >> (interpC[0] - 1));
            }
            // 二次 B 样条，在整数区间中点插值
            if (Ord == 2 && type == BCType::notAknot)
            {
                // 注意一共有 N 个插值点，但两端都是整数点，需要排除
                // 插值点的个数
                N = N - 1;
                vector<double> b(N - 1);
                vector<vector<double>> M(N - 1, vector<double>(N - 1));

                //系数多项式
                for (int i = 1; i < N - 2; i++)
                {
                    M[i][i - 1] = M[i][i + 1] = 1;
                    M[i][i] = 6;
                }
                M[0][0] = M[N - 2][N - 2] = 5;
                M[0][1] = M[N - 2][N - 3] = 1;

                //目标向量
                for (int i = 1; i < N - 2; i++)
                {
                    b[i] = 8 * interpC(interpC[i + 1]);
                }
                b[0] = 8 * interpC(interpC[1]) - 2 * interpC(interpC[0]);
                b[N - 2] = 8 * interpC(interpC[N - 1]) - 2 * interpC(interpC[N]);

                //获得样条系数
                vector<double> a(N + 1);
                vector<double> tmp_a = Solver::solve(M, b);
                for (int i = 0; i < N - 1; i++)
                {
                    a[i + 1] = tmp_a[i];
                }
                a[0] = 2 * interpC(interpC[0]) - a[1];
                a[N] = 2 * interpC(interpC[N]) - a[N - 1];

                //累计样条函数
                res = BSpline<2>(0) * a[0];
                for (int i = 1; i < N + 1; i++)
                {
                    res = res + BSpline<2>(i) * a[i];
                }
                //根据左端点的位置平移
                res = (res >> (interpC[0] - 1));
            }
            break;
        }
        }

        return res;
    }

    // 二维曲线拟合
    template <int Ord>
    Spline<2, Ord, SplineType::ppForm> fitCurve(const vector<Vec<double, 2>> &points, BCType type)
    {
        // 通过累计弦长转化为一维情形
        vector<double> lengths(points.size());

        // 计算累计弦长
        lengths[0] = 0;
        for (int i = 1; i < points.size(); i++)
        {
            double dx = points[i][0] - points[i - 1][0];
            double dy = points[i][1] - points[i - 1][1];
            lengths[i] = lengths[i - 1] + sqrt(dx * dx + dy * dy);
        }

        //分别计算插值多项式
        InterpC interpC_x;
        InterpC interpC_y;
        interpC_x.points = lengths;
        interpC_y.points = lengths;

        for (int i = 0; i < lengths.size(); i++)
        {
            interpC_x.interpSites[lengths[i]] = {points[i][0]};
            interpC_y.interpSites[lengths[i]] = {points[i][1]};
        }

        Spline<1, Ord, SplineType::ppForm> res_spline1 = interpolate<Ord, SplineType::ppForm>(interpC_x, type);
        Spline<1, Ord, SplineType::ppForm> res_spline2 = interpolate<Ord, SplineType::ppForm>(interpC_y, type);

        // 赋值给结果多项式
        Spline<2, Ord, SplineType::ppForm> res_spline(lengths);
        for (int i = 0; i < res_spline.size(); i++)
        {
            // 合并多项式为二维多项式
            res_spline[i] = overlying<Ord, 1, 1>(res_spline1.spline_poly[i], res_spline2.spline_poly[i]);
        }
        return res_spline;
    }
};

// 为了防止宏的重复定义
#undef MAX
#undef MIN